# A2 Density Plot

# Load packages 
library(readstata13)
library(ggplot2)
library(dplyr)

# Read Replication data 
data <- read.dta13("~/Dropbox/Replication_MVC/Datasets/datasets_analysis/panel_excombatientes.dta",
                   convert.factors = TRUE, generate.factors = FALSE,
                   encoding = "UTF-8", fromEncoding = NULL, convert.underscore = FALSE,
                   missing.type = FALSE, convert.dates = TRUE, replace.strl = TRUE,
                   add.rownames = FALSE, nonint.factors = TRUE, select.rows = NULL)

# Captures by mun per 25 ex-combatants
result <- data %>%
    arrange(origmun, year) %>%
    group_by(origmun, year) %>%
    summarise(count_id = n(),
              sum_diffcaptures = sum(diffcaptures))

result$arrestsper25excomb <- result$sum_diffcaptures*25/result$count_id

# Separate for each year
attach(result)
result_2014 <- result[ which(year=='2014'),]
detach(result)

attach(result)
result_2015 <- result[ which(year=='2015'),]
detach(result)

attach(result)
result_2016 <- result[ which(year=='2016'),]
detach(result)

density_14 <- ggplot(result_2014, aes(x=arrestsper25excomb)) + 
    geom_density(fill="blue", alpha=.2) +
    geom_line(stat="density")+
    xlim(1, 28)+
    xlab("")+
    ylab("") +
    ggtitle("(2014)")

density_15 <- ggplot(result_2015, aes(x=arrestsper25excomb)) + 
    geom_density(fill="blue", alpha=.2) +
    geom_line(stat="density")+
    xlim(1, 28)+
    xlab("Arrests per 25 Ex-Combatants")+
    ylab("Density") +
    ggtitle("(2015)")

density_16 <- ggplot(result_2016, aes(x=arrestsper25excomb)) + 
    geom_density(fill="blue", alpha=.2) +
    geom_line(stat="density")+
    xlim(1, 28)+
    xlab("")+
    ylab("") +
    ggtitle("(2016)")

multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
    library(grid)
    
    # Make a list from the ... arguments and plotlist
    plots <- c(list(...), plotlist)
    numPlots = length(plots)
    
    # If layout is NULL, then use 'cols' to determine layout
    if (is.null(layout)) {
        # Make the panel
        # ncol: Number of columns of plots
        # nrow: Number of rows needed, calculated from # of cols
        layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                         ncol = cols, nrow = ceiling(numPlots/cols))
    }
    
    if (numPlots==1) {
        print(plots[[1]])
        
    } else {
        # Set up the page
        grid.newpage()
        pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
        
        # Make each plot, in the correct location
        for (i in 1:numPlots) {
            # Get the i,j matrix positions of the regions that contain this subplot
            matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
            
            print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                            layout.pos.col = matchidx$col))
        }
    }
}

Density_13_16 <- multiplot(density_14, density_15, density_16, cols=3)
